function [IdxTree, ErrTree, FlowRec, History, IdxDict] = swincresym(A, c)
% Smallwhite Symmetric NMF Package (swsymnmf)
% by Sun Sibai (Smallwhite <niasw@pku.edu.cn>)
%
% This is the incremental function of Symmetric NMF. The algorithm is ANLS.
%
% Variables:
%   >> A: sparse symmetric matrix to be factorized (Link Matrix)
%   >> c: number of dimensions of factors (Cluster Number)
%   << IdxTree: cluster result of each step (Cluster Result)
%   << ErrTree: residual error of each step (to be minimized)
%   << FlowRec: maxflow record of each step (to be minimized)
%   << History: cluster result of each step (Historical Clusters in Json)
%   << IdxDict: cluster result of each step (Cluster Result in Json Dict)
%
% References:
%   * symnmf_anls [http://www.cc.gatech.edu/~dkuang3/]
%   > by `Da Kuang`, `Chris Ding`, `Haesun Park`,
%   > Symmetric Nonnegative Matrix Factorization for Graph Clustering,
%   > The 12th SIAM International Conference on Data Mining (SDM '12), pp. 106-117.
%
% Basic Mechanism:
%   W = (As'As+Ad Ad');
%   min_{H,W} er(H) = ||M - WH'||_F^2 + a * ||W-H||_F^2
%   H_{j,k}>=0, a = max(M_{j,k})
%
%
 c=floor(c);
 if (c<=1)
  error('cluster number c should > 1 and be integer.');
 end
 if (size(A,1)~=size(A,2))
  error('adjacent matrix A should be a square matrix.');
 end
 if (c>size(A,1))
  error('cluster number c should < dimension of adjacent matrix. (paper number > cluster number)');
 end

 IdxTree = zeros(size(A,1),2);
 ErrTree = zeros(c-1,2);
 FlowRec = zeros(c-1,2);
 History = {};
 IdxDict = {};
 IdxTree(:,1) = 1:size(A,1);
 Cluster = unique(IdxTree(:,2));
 History.c1 = Cluster;

 Mt = A*A'+A'*A;

 for itstep = 1:(c-1)
  flowMax = 0; % which seperating maximizing flows
  clstDiv = 0; % which cluster should be seperated
  IdxTreeNew = IdxTree;
  for itclst = 1:size(Cluster,1)
   disp(sprintf('trying: step(%d/%d),cluster(%d/%d)',itstep,c-1,itclst,size(Cluster,1)));
   subidx=find(~(IdxTree(:,2)-Cluster(itclst))); % the subset
   As = A(subidx,:);
   Ar = A(:,subidx);
   %M=As*As'+Ar'*Ar-(Al*Al'+Al'*Al)/2; % similarity
   M=Mt(subidx,subidx);
   Al=A(subidx,subidx);
   M=M-(Al*Al'+Al'*Al)/2;
   W=swnormal(M);
   [H,iter,err]=symnmf_anls(W,2,struct('debug',1)); % ANLS symnmf
   H=full(H);
   [clstmp,maptmp]=sort(sum(H)); % sort by cluster size
   Hs=zeros(size(H));
   for it=1:size(maptmp,2)
    Hs(:,it)=H(:,maptmp(it));
   end
   [Max,Idx]=max(abs(Hs'));
   Idx=Idx+2*(Max==0); % 3 for unrelated cluster
   Idx=Idx';
   IdxTreeTry=IdxTree(:,:);
   IdxTreeTry(subidx,2)=IdxTreeTry(subidx,2)*4+Idx; % quaternary ID
   flowMaxTry=sum(sum(swflow(A,IdxTreeTry).^2)); % just sum up (may replaced by fnorm)
   if (flowMaxTry>flowMax)
    flowMax=flowMaxTry;
    clstDiv=Cluster(itclst);
    ErrTree(itstep,:)=[clstDiv,err];
    FlowRec(itstep,:)=[clstDiv,flowMax];
    IdxTreeNew=IdxTreeTry;
   end
  end
  IdxTree=IdxTreeNew;
  Cluster = unique(IdxTree(:,2));
  History = setfield(History,strcat('c',sprintf('%d',itstep+1)),Cluster);
 end
 disp('Creating Dictionary ...');
 for it=1:size(IdxTree,1)
  IdxDict = setfield(IdxDict,strcat('p',sprintf('%d',IdxTree(it,1))),IdxTree(it,2));
 end
 disp('Done.');
end
